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Abstract 

We propose a new model for the description of complex granular particles 
and their interaction in molecular dynamics simulations of granular material 
in two dimensions. The grains are composed of triangles which are connected 
by deformable beams. Particles are allowed to be convex or concave. We 
present first results of simulations using this particle model. 

I. INTRODUCTION 

Molecular dynamics has been proven to be a well suited method to investigate the dy- 
namic and static behavior of granular matter. The concept of molecular dynamics was 
initially used to simulate the dynamics of atoms and molecules, i.e. of particles without 
inner degrees of freedom. Alder and Wainwright who belong to the inventors of this me- 
thod investigated already 1957 numerically the low density approximation of the entropy 
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S = —N ks (In f(p,t)) (/ is the one-particle probability density in momentum space) in a 
hard sphere system consisting of N = 100 particles [1]. There are applications of molecular 
dynamics in many fields of interest and molecular dynamics became one of the standard 
methods in computational physics. Recent molecular dynamics simulations are of very high 
algorithmic complexity, and systems consisting of up to 10 9 particles (with simple interaction 
forces) have been simulated (see e.g. [2-4]). A very interesting overview on the historical 
development is given in the first part of Hoover's book [5] . 

In molecular dynamics of granular particles, sometimes called "granular dynamics", 
each of the "molecules" is a macroscopic body typically with a diameter of the order 
D 0.1... 1 mm which has its own thermodynamic properties, i.e. it has internal de- 
grees of freedom and hence it can dissipate mechanical energy when it is subjected to outer 
forces. In some special cases as in simulations of planetary rings [6,7] the diameter can also 
be of the order D « 1 mm ...10m. In the most simple case one can express the inelastic 
nature of the collisions of the grains by introducing restitution coefficients which describe 
the relative velocities in normal and tangential direction vff and of two granular particles 
i and j after a collision as functions of the relative velocities before the collision and 
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Here the relative velocity is 
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or, including the particle rotation with angular velocities Ui and ujj , by 
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Hence the relative velocities at the point of contact are 
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In general the normal and tangential restitution coefficients and are themselves 
functions of the relative impact velocities e N = e N (v^ , V^j) and e T = e T (v^ , V?^. For 
the case that the particles are rough homogeneous spheres the coefficients and can 
be calculated analytically [8] solving the viscoelastic equations of the particles. Using this 
model one neglects the complicated interaction between particles during the collision. The 
approximation is justified if collisions are relatively rare events, i.e. if the mean free path 
between collisions is long and one can restrict the calculation to two-particle-interactions 
only. If a system is in this regime one calls it a "granular gas". Event driven algorithms 
have been applied in many papers, recently e.g. in [9,10] for an investigation of inelastic 
clustering of granular particles (using constant restitution coefficients) and by Luding et 
al. to simulate vertically shaken material [11,12]. For sophisticated implementations of such 
algorithms see [13,14]. 

In many cases of interest, however, the particles collide very often or they touch each 
other permanently, i.e. the system reveals static properties. Examples for such systems are 
flows through pipes and hoppers, ball mills, shaken materials, sand heaps and many others. 
In the model proposed by Cundall and Strack [15] and Haff and Werner [16] the particles 
feel restoring forces in normal and tangential direction, and they lose kinetic energy during 
the collisions. If the particles i and j at the positions fi and fj with radiuses Ri and Rj 
touch each other they feel the force 
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Y is the Young modulus, 7^ and 7^ are the damping coefficients in normal and tangential 
direction and \i is the Coulomb friction coefficient. Eq. (6d) describes the relative velocity 
of the surfaces of the particles at the point of contact and eq. (6e) gives the effective mass. 
Eq. (6c) takes the Coulomb friction law into account, saying that two particles slide on top 
of each other if the shear force overcomes \x times the normal force. Eq. (6b) comes from the 
Hertz law [17] for the force between two rigid spheres which are in contact with each other. 

This model has been successfully applied to simulate the behavior of dry granular ma- 
terial in many works describing various physical phenomena. It developed to the standard 
method for calculations of granular dynamics. Some of these problems are size segregation 
and convection in vibrated material in two and three dimensions (e.g. [18-21]), the flow in 
hoppers (e.g. [22-24]) and pipes (e.g. [25,26]), the flow in rotating cylinders (e.g. [27-29]), 
the motion of granular material on an inclined surface [30,31], sound propagation in granular 
material [32], the onset of turbulence [33,34] and many others. Many experimental results 
of many authors could be reproduced numerically using this type of molecular dynamics 
simulations. There have been developed very efficient algorithms for the case of short range 
interaction with large force gradients as it is typical for granular materials, e.g. [35-37]. The 
model by Cundall and Strack yields reliable results in more dynamic systems, i.e. when the 
static friction of the particles does not play a major role. When the static properties of 
the material govern the behavior, however, the model might fail. In [38] is shown that one 
cannot build up a stable sand heap with a finite inclination of such particles but the heap 
dissolves under its own mass and the angle becomes smaller with rising number of particles. 
Lee [39] simulated static friction effects using spheres which are connected by springs. When 
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two particles touch each other for a certain time a spring "grows" between these particles 
and keeps the particles together when forces are applied which would separate them. When 
the separating forces become too large the springs break and the particles move freely. Using 
this model Lee reproduced the finite angle of repose in a sand heap, however, the growing 
and breaking springs have no direct equivalent in nature, and hence the model seems to be 
quite artificial. 

Another simulation method which corresponds to the previous one for the case of very 
large damping has been introduced by Visscher and Bolsterli [40]. In each iteration step all 
of the grains are moved according to the motion of a vessel containing the granular material. 
Then the particles are released one by one, beginning with the lowest one (with respect to 
the direction of gravity). Each particle moves until it reaches the next local minimum of 
the potential energy when it hits either the wall or other previously released grains. When 
a particle reached its local minimum it stays in that position, i.e. a falling particle cannot 
cause the motion of another previously released particle. This behavior means that inertia 
of the moving particles is neglected as soon as they hit another particle or the wall, and 
hence it corresponds to the case of very large damping. The advantage of this algorithm 
is its numerical simplicity, i.e. it is possible to investigate much larger systems than using 
molecular dynamics. Recently it has been applied in simulations of several problems as size 
segregation [41] and the motion of particles in a rotating cylinder [42,43]. Obviously the 
simplification of infinite damping is not valid in each case, and hence one has to carefully 
investigate whether it is justified to apply this algorithm (see e.g. [45]). 

There are some more simulation methods for the dynamics of granular material which we 
want to mention only: Peng and Herrmann applied a Lattice Gas Automaton to investigate 
density waves in a pipe [46]. Caram and Hong investigated granular flows using a Random 
Walk technique [47] . 

There are some models to simulate complex shaped grains which are composed of spheres. 
In the model by Gallas and Sokolowski [48] the grains consist of two spheres connected with 
each other by a stiff bar. Walton and Braun applied more complicated particles consisting of 
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four or eight spheres rigidly connected with each other [49] . Using this model they examined 
the transition from stationary to sliding flow and the transition from sliding to raining flow 
in a rotating drum. Poschel and Buchholtz [50,38] describe grains built up of five spheres 
where one of them is located in the center of the grain and four identical spheres are at the 
corners of a square. Each pair of neighboring spheres is connected by a damped spring. The 
latter model was applied to the rotating cylinder [29] and it was shown that the simulation 
results agree much better with the experiment than equivalent calculations with spheres. 
Especially it was shown that one can reproduce stick-slip motion and avalanches which was 
not possible in simulations using spheres. The inclination of the surface of the material in 
the rotating cylinder and the dependency of the inclination on the angular velocity, however, 
did not agree well with the experiment. 

Mustoe and DePorter [51] proposed a particle model where the boundary geometry of 
the particles is defined (in local coordinates) by 

ni 

-1 = 0. (7) 

Varying rii between 2 and oo the shape of the ith particles varies continuously from elliptical 
to rectangular. To detect contacts of pairs of particles described by eq. (7) one has to 
solve numerically equations of high order for each pair of possibly touching particles in 
each iteration step. This requires an iteration technique which converges slowly for higher 
exponents rii and which makes hence the algorithm numerically complicated. Hogue and 
Newland [52] investigated the flow of granular material on an inclined chute and through a 
hopper using a model where the boundaries of the convex particles are given by polygons 
with up to 24 vertices each. To detect whether two particles touch each other one has 
to calculate the intersections between each pair of vertices. During collisions energy is 
dissipated according to Stronge's energy dissipation hypothesis [53] in normal direction, 
whilst Coulomb's friction law models the energy losses also in tangential direction. 

Using grains which consist of interconnected spheres or of particles described by eq. (7) 
it is not possible to simulate particles with sharply formed corners. For some effects it seems 
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to be essential to simulate such particles to reproduce the experimental observed effects. 
This point is discussed in detail in [38]. Tillemans and Herrmann [54,55] proposed a two 
dimensional particle model where the particles are convex polygons. When two particles i 
and j touch each other, i.e. when there is an overlap, there acts the force 



A is the compression area of the particles (overlap area), Y is the Young module of the 
material. The effective mass ra^ is given by eq. (6e), the relative particle velocity is v rel = 
fi — fj, where fi points to the center of mass of the particle i, and L c is the characteristic 
size of the particle. The unit vectors are in the direction of the line that connects the 
intersection points of the overlapping particles i and j (e T ) and perpendicular to this line 
(e N ). Only convex particles are allowed. This model was used e.g. in simulations of shear 
cells, earth quakes and flow through hoppers [55]. The results differ significantly from similar 
calculations using spheres. For the case of the hopper simulation they found clogging and 
arching which could not be found in simulations with spheres. A similar model was used 
earlier by Handley [56] who investigated the fracture behavior of brittle granular material. 
Potapov et al. [57,58] investigated a similar model for solid fracture. Initially they subdivide 
a macroscopic two dimensional body into many small equilateral triangles (elements). The 
forces in the body are resolved as forces at inter-element contacts of two different types 
which they call "glued" and "collisional" . Glued contacts are the joints between elements 
interior to a solid body which can support stresses. When the tensile stresses exceed a certain 
value, cracks form and the glued bond breaks. Collisional contacts are contacts between the 
surface elements of the body (which eventually collides with other bodies or walls), and 
contacts between triangles in the bulk of the body where glued contacts have been broken. 
In the three latter models [54-58] the convex grains have been Voronoi polyhedra in two 
dimensions [59]. 
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In this paper we present a new particle model for molecular dynamics of granular material 
in two dimensions where the particles are simulated using a similar model which has some 
advantages. In the following section we describe the model, the forces acting between the 
grains are derived in sections III and IV. Some aspects of the implementation and the 
performance of the algorithm are discussed in section V and some sample results based on 
this particle interaction model will be briefly discussed in section VI. 

II. THE MODEL 

In our model the grains are composed of an arbitrary number of ideal elastic triangles 
which are connected by beams (see fig. 1). A beam in our sense is an deformable damped 
bar which is subjected to forces in the direction of its axis and transverse to its axis, i.e. to 
normal and shear forces, and to moments acting on its ends. Triangles which belong to the 
same grain do not interact with each other. When two triangles of different grains collide, 
i.e. if there is an overlap between both, they feel a restoring force. Hence the force acting 
upon the triangle i belonging to the grain j (j G {1,2,..., N}) is 

JV n t (Z) n{(i) 

f/= £ £r2+£*£ + *(^)- (9) 

l=l,l^j k=l fc=l 

The first sum in the first term in eq. (9) runs over all grains / except the jth, the second sum 
runs over all triangles k G {1,2,..., n t (l)}, the lih grain consists of. ff k is the force which 
acts between the triangles i and k which belong to different grains j and / (see section III). 
The sum in the second term runs over all beams connecting the triangle i of the grain j 
with other triangles of the same grain. A/ fc is the force induced by the beam k acting on the 
ith triangle of the grain j. They originate from the distortion of the beam that connects 
the triangle i with another one of the same grain. The acting forces and momenta will 
be discussed below in section IV. A similar model for the beam forces has been introduced 

~~ * '3 

earlier by Herrmann et al. for the fracture of disordered lattices [60]. The last term , r i ) 
describes the action of an external force, as e.g. gravity, on the triangles. 
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During its deformation a beam dissipates energy similar to a linearly damped spring, 
i.e. proportional to its deformation rate (see section IV). To avoid time consuming calcu- 
lations due to Steiner's law each beam is fixed at both ends in the center of mass points 
of the triangles it connects. When the beam is not deformed by eventually applied forces 
(the beam is at rest) it is perpendicular to the neighboring edges of the triangles which it 
connects. Fig. 1 shows some examples of grains. The only restriction concerning the num- 
ber, the shape and the position of the triangles is caused by the condition that the beams 
are fixed in the center of mass. 

The next section III gives a detailed description of the force that acts when triangles 
collide, in section IV the forces due to deformed beams are derived. 

III. THE INTERACTION OF THE TRIANGLES 

In this section we describe the calculation of the term Tf k in eq. (9) originating from the 
compression of two triangles i and k which belong to the grains j and /, respectively. When 
the triangles are compressed there exists a virtual "overlap area" which leads to an elastic 
restoring force T ik . For simplicity here and in the following we drop the upper indexes of 
the variables. Areas of triangles XYZ and quadrangles WXYZ we denote by A( XYZ) and 
0(WXYZ). Fig. 2 shows the most frequently found type of a collision between the triangles 
AiBiCi and A k B k C k , i.e. a corner of one triangle deforms a side of another one. The absolute 
value of the restoring force is given by the shadowed area A (SiS 2 Ak) of the triangle 
SiS 2 A k times the Young module Y. Its direction is perpendicular to the intersection line 
SiS 2 (Poisson hypothesis, see e.g. [23]). Hence resulting momenta acting upon the triangles 
AiBiCi and A^B^C^ read 

M (A t B t Q) = HGi x f lk (10a) 
M (A k B k C k ) = HG k x f fci = -H~G k x f ik . (10b) 

The vector HGi points from the middle point of the line 5152 to the center of mass of the 
triangle AiBiCi. 
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Although the case shown in fig. 2 is the mostly occurring type of collisions there are 
several other types which treatment has to be discussed below. 

Classifying the possible interactions of the triangles we define five different types of col- 
lisions which are slightly different from those in [57]. The detection of overlapping polygons 
is a very important problem in computer graphics too, where one frequently has to decide 
whether two objects cover each other, and which of them is in the foreground or in the back- 
ground, respectively. Therefore we found it helpful and inspiring first to have a look into 
advanced methods in computer graphics, e.g. [61,62]. In all figures 2-7 the overlapping area 
is drawn extremely exaggerated. Typically the overlapping area of two interacting triangles 
does not exceed a few tenth of a percent of the area of the triangles. 

1. For the first type (fig. 3) there are two intersection points Si and S 2 lying at the same 
edge of one of the interacting triangles. One of the situations drawn in fig. 3 corres- 
ponds to the standard collision (fig. 2). We calculate the area A (SiS 2 Ai) (shadowed 
area) given by the intersection points Si and S 2 and the included point A{. Because 
the overlapping area is definitely smaller than half the area of the triangles the force 
is proportional to the minimum between A {S\S 2 Ai) and A (AiBiC'i) — A (SiS 2 Ai) 



The force acts perpendicular to the line between the intersection points S1S2. 

2. The second type (fig. 4) is characterised by two intersection points lying on diffe- 
rent edges for both triangles. We calculate the areas A (SiS 2 Ai) (dark gray) and 
A (SiS 2 Ak) (light gray). The force is proportional to the overlap 



It acts perpendicular to the line between the intersection points as in the previous 




T tk = Y ■ min {A (S^A) , A (A t B t C t ) - A (S^A,)} . 



(11) 




T lk = Y- [min {A (S^A,) , A (A t B t Q) - A (S^A,)} + 



min{A (S^A,,) , A (A k B k C k ) - A (S^A,)}] 



(12) 



case. 
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3. The third type (fig. 5) is characterised by four intersection points, each pair of two 
points lies on one edge of both triangles. Therefore the third edge of both interacting 
triangles does not intersect the edges of the other triangle. There are two pairs of 
forces proportional to the quadrangular area □ (S1S2S3S4) of the overlap 



ik 



T ik = \ ■ Y • □ (SiS 2 S 3 S 4 ) , (13) 



one acting perpendicular to the line SiS 2 and the other perpendicular to the line S1S4. 
The pre-factor 0.5 is due to two pairs of acting forces instead of one for the first two 
collision types. 

4. The fourth type (fig. 6) is characterised by four intersection points too, but here one of 
the triangles has intersection points at all its edges. Because interactions of this type 
occur extremely seldomly we did not implement an exact calculation for this case, but 
we calculate two interactions of type 1 with the intersection points Si-S 2 or 5 3 -5 4 
respectively, instead of solving the correct problem. In real simulation the fraction of 
this type occurs approximately 1CT 5 of the number of collisions. Nevertheless one has 
to deal with these extremely seldom events since otherwise they might cause problems 
from which the system cannot recover. If one does not care about these events usually 
the system gains in a few time steps after the event occurs a huge amount of kinetic 
energy, i.e. the system explodes. 

5. The fifth and last type of interaction (fig. 7) is characterised by six intersection points 
Si - Sq. Equivalent to the forth type we do not calculate the exact interaction but 
substitute this calculation by solving three interactions of the first type with pairs of 
interaction points Si-S 2 , S 3 -S 4 and S 5 -S 6 . 

According to the forces we derive forces in parallel to the axes of the Cartesian coordinate 

— * 

system and moments M acting with respect to the center of mass points Gi and Gk as 
described in equations (10). 
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IV. STRESSES IN BEAMS 



When torques and forces are applied to a beam the beam deforms. Since all deforma- 
tions are assumed to be small compared with the size of the beam one can superpose the 
deformations originating from different applying forces and torques [63,64]. Fig. 8 shows an 
(infinitesimally) deformed beam with radius of curvature p. We find the approximation 

f d<d\ 1 dx , H A . 

tan — = • dx or p = — — 14 

V 2 J 2 ■ p ' dG K J 

with tan(dG) = dO. 

The dashed line ss is the neutral surface length of which does not change during the 
deformation. All other fibres below and above the neutral surface are in tension or in 
compression, respectively. The length of a fibre in distance y from the neutral surface is 
(p + y) dO = (1 + y/p) dx and hence the strain reads e = y/p. With Hooke's law a = E ■ e 
one finds a = Ey/p. Finally the resulting momentum M is 

EI 



/ oydA = — (15) 

J A P 



M 

P 

where the moment of inertia I = J y 2 dA depends on the geometrical shape of the cross 
section of the beam. In the present paper we assume that the beams have quadratic cross 
section of width b and hence / = 6 4 /12. When we express the radius of curvature p by 

1 d 2 v 

r I? • (16) 

where v is the deflection of the beam from its initial position (with the approximation 
tanG Rj G). From eqs. (15) and (16) one finds the basic equation for the deflection of a 
beam: 

d 2 v „ M 

l^ = V = + EI (17) 

This equation has to be solved to find the forces and moments which act when a beam 
is deformed. The integration constants are used to satisfy the boundary conditions. 
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There are three different deformation modes for beams: elongation, shearing and ben- 
ding. In the following three subsections IV A, IV B and IV C we will discuss the deformation 
rules in detail. 



A. Elongation of beams 

Since the exact notation of the acting forces and moments as vector functions of the 
coordinates of the triangles the beam connects, is not very instructive but confusing due to 
the length of the formulae expressions we discuss for simplicity of notation the deformation 
of the beams here and in the following in local coordinates. Instead of providing the exact 
vector notation we rather discuss absolute values. In all cases the directions in which the 
forces and moments act are very clear, moreover they are given for each case explicitly in the 
figures 9-11 and in figure 14. For the implementation of the algorithm, however, one should 
note that it is necessary to transform the given expressions for the forces and momenta into 
the coordinate system of the triangles. These transformations require a considerable part of 
the computation time. 

Fig. 9 displays the shape and the location of a beam of length L at rest, i.e. when no 
deforming forces Fa, F b or moments Ma, M b at the ends A and B of the beam apply. The 
ir-axis is drawn with a thin solid line, i.e. at rest the beam lies on the ir-axis. 

For the elongation deformation we find according to the linear Hooke law the restoring 
force 

F el = |AL| • E. (18) 



B. Shearing 

We want to calculate the moments and forces M A , M B , F A and F B of a sheared beam as 
drawn in Fig. 10. For a beam of length L/2 which is fixed at one end (at x = 0) and where 
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acts the shear force F (in negative ^-direction) at the other end (at x = L/2) we find for 
the local moment M(x) = F — x^j and hence eq. (17) reads 

v" = (- - x] (19) 



EI \2 

with the boundary conditions v(0) = and v'(0) = 0. We find for the vertical deviation 
A/2 at the point x = L/2 

FL 3 

A/2 = — 20 

' 2AEI v 1 

force which acts at the free end of the beam. 

One can assume that the sheared beam in Fig. 10 is built up of two of such beams of 
length L/2 as described above. Finally one finds 

FL 3 , N 

A = 21 

12 EI y ' 

and hence 

1 2FT 

F A = -F B = F = — A (22a) 

M A = M B = F ■ | = ^A . (22b) 



C. Bending 

When a beam undergoes bending deformation (fig. 11), the resulting forces and moments 
can be found by superposing the forces and moments that would act if the beam would be 
bent at each side separately. Therefore we can restrict ourself on calculating the forces and 
moments according to the single sided deformation drawn in fig. 12. 

The deformation in fig. 12, however, can be understood as a deformation of a beam 
with free ends according to a single sided moment M (fig. 13) superposed by a moment M b 
(fig. 12) that assures that the angle Q B (fig- 12) results to zero. 
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First we consider one sided bending of the beam (fig. 12). For a beam as drawn in fig. 13 
where the moment M acts at one of its both sides one finds easily 

ei = (23b) 

B 6EI v ; 



Since the beam in Fig. 12 has the angle Q B = we find the moments acting at the ends 
of the beam in fig. 12 by superposing the moment M B , which causes a virtual angle Q* B , 
fulfilling the condition G^* + 0* B = 0: 

m b = 5f£ ( -ey = f . (24) 



The moment M B causes the angle 



Ua ~ 6EI ~ 12EI {Zb) 
and the resulting angle G^ (fig. 12) is 

e A = e\ + &X = ^. (26) 

Hence, finally we find for the different moments and forces for the bending deformation 
(fig. 12) 

4EI 

M A = - — e A (27a) 
2EI 

M B = — j-@ A (27b) 

and with 

M A = -F B L - M B (28) 

F B = -F A = -jy®a ■ (29) 

(Note that the moments M A and M B in eq. (28) have to be added as if they would apply to 
the same end of the beam.) 
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If the beam is bent at both ends by angles Q A and Q B (Fig. 11) the resulting moments 
M A and M B and forces F A and F B can be calculated by a superposition of the forces and 
moments according to two independent bending deformations of the type discussed in the 
previous case. 

The angle 0^ generates the moments M* A and M B according to eqs. (27) 

4£J 



-0 A (30a) 



L 
2EI 

M* B = — S A (30b) 



and the angle Q B causes the moments 



2EI 

MX = — —Q B (31a) 
4EI 

M B * = — e B . (31b) 



Hence we find the resulting moments 



/4EI 2EI \ 
M A = M* A + M* A * = - [—Qa + —®b) (32a) 

/2EI 4EI \ 
M B = M* B + M* B * = - l—e A + —e B , (32b) 



V L 

and with M A + M B = —F B L the resulting forces 

QE I 

f b = -f a = — (e A + e B ) . (33) 

Comparing the results from sections IV B and IV C one can show that each deformation 
of a beam can be expressed by bending and elongation, i.e. the shear deformation need not 
to be considered. The proof is given in Appendix A. 

The deformation of the beams is damped proportionally to the deformation rates 



Mf = -21 e A , 


(34a) 




(34b) 


F A ^ = —7 (v A — v B ) ■ AB, 


(34c) 


F^ d) =1 (v A -v B )-AB, 


(34d) 
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where 7 is the damping coefficient of the beam material. As proven in the appendix the 
deformation of the beam is in linear approximation completely determined by the angles 0^ 
and Qb and the length L. Hence there is no damping force acting in shear direction. 

V. IMPLEMENTATION OF THE ALGORITHM 

We implemented a molecular dynamics algorithm using the particle model described 
in the previous sections in FORTRAN. A Gear predictor corrector scheme of fourth or- 
der [65,66] was applied to integrate Newton's equations of motion. We have to proceed both 
the forces caused by the elastic deformation of the triangles as well as the forces induced 
by the beams for each iteration step (eq. 9). Since every particle can interact with each 
other one causing a high algorithmic complexity, we applied a Verlet neighborhood list [67] 
to decrease the amount of computer time. As discussed below (see table I) the calculation 
of the triangle intersections is the most time intensive part of the calculation. Hence we 
applied an advanced version of the neighborhood list method to keep the fraction of the 
computation of the triangle-interactions as small as possible. 

Our list method acts in two steps. In the first step we prepare neighborhood lists for 
each triangle to reduce the number of possible interactions which has to checked. Because 
we have to investigate all possible triangle interactions the time for this step rises as the 
square of the total number of triangles T ~ (N ■ n t ) 2 . (For simplicity we assumed for 
this estimate that all N grains consist of n t triangles each. Hence there are (N ■ n t ) ■ 
((N — 1) • n t ) possible triangle-interactions.) In the second step we check for each entry of 
the constructed neighborhood list whether the corresponding triangles interact indeed or 
not. If they interact, i.e. if there is an overlap area, the type of interaction is determined 
according to the classification scheme in section III. 

Finally we construct three lists for the interactions of the first three types of the clas- 
sification scheme. The interactions of the forth and fifth type are replaced by interactions 
of the first type as described above. The computer time for this step rises linearly with the 
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number of particles N and with the number of triangles n t belonging to each particle. Dif- 
ferent from usual Verlet tables these lists contain only triangles which definitely do interact, 
and hence we call the most time consuming subroutines for the calculation of the forces only 
for triangles which do interact in one of the five manners described above (section III). 

The calculation of the forces induced by the beam deformations follows straight forward 
the formulas described in section IV. 

In table I we present a detailed performance analysis of the algorithm. The numbers in 
the table refer to simulations on a DEC 3000/700 workstation. As visible from the data most 
of the time (approximately 70%) is necessary to construct the interaction list, because we 
check here every pair of neighbored triangles for intersections. The construction of this lists 
give the possibility to classify the interactions and therefore to use for each kind of interaction 
its own optimised algorithm. This prevents us from calculating useless intersection points, 
areas and other data. 

The computer time needed for the predictor-corrector integration is very low (3.7%). 
Thus it seems to be not useful to optimise the integration procedure, or to decrease the 
accuracy of the integration to safe computer time. The calculation of the beam forces is not 
time expensive as well. 

From the data in table I it is obvious that triangle interactions of the first type domi- 
nate. Approximately 98% of the occurring interactions are of this type. Interactions of the 
fourth and fifth type, where our program yields approximative results only, occur extremely 
seldomly. We found them only a few times during all our simulations. 

VI. FIRST RESULTS 

A. The collision of a two particles 

To demonstrate the non-trivial behavior of colliding non-spherical particles and to check 
the correctness of our implemented model we want to present the results of an experiment 
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where a quadratic grain collides with another resting one in the absence of gravity. The 
grains consist of four equal triangles each, forming a simple quadratic grain as shown in 
fig. 1. The parameters of the materials are: Y = 2 • 10 7 g/ (cm sec 2 ), E = 1 • 10 5 g /(sec 2 ), 
1 = 1- 10~ 4 cm 3 and 7 = 9 gj sec. The value of the Young-module Y agrees with most of 
the simulations of spheres which can be found in the literature. 

Figs. 15 (left figures) show the time evolution of the particles as a sequence of stro- 
boscopic snapshots for different relative velocities of the particles a) v re i = 10 cm/ sec, 
b) v re i = 50 cm/sec and c) v re i = 100 cm/sec. One of the triangles of each grain has 
been filled to visualise the rotation of the grains. Initially the velocity of the left particle 
is v = (v re i,0) cm/sec and the right particle rests. For low velocities v < v tc ~ IQcm/sec 
the grains collide twice within a very short time. Therefore the traces of the particles 
that undergo soft collisions with velocity values v < v tc differ qualitatively from the 
traces of particles colliding with higher relative velocity. Since it is hard to display the 
complicated motion of the grains in detail we refer to our World Wide Web-site URL 
http://summa.physik.hu-berlin.de:80/~thorsten/MDASGP.html [68] where one can find 
animated sequences of various problems mentioned in the current paper. The right hand 
side in fig. 15 shows the rotational and kinetic energies and the total energy as a function 
of time. The initial kinetic energy of the translational motion of the left particle results 
in kinetic energy of both particles due to translational motion and rotation. A part of the 
mechanical energy is lost due to dissipation in the beams. The relative amount of rotatio- 
nal energy depends strongly on the impact velocity as expected from the discussion above. 
The dependence of the extinction of the rotational degree of freedom on the impact rate is 
demonstrated in fig. 16 too. This figure shows the time evolution of the angular velocity 



of the grains to = — 



nt 



E (Vi ~ V g ) X S, 



=1 



over time for different values of the impact velocity. 



v g = — Y, Vi is the velocity of the grain and Si is the vector pointing from the center of mass 



i=l 



point of the grain to the center of mass point of the triangle i. n t is the number of triangles 
the grain consists of. For high velocity one observes an abrupt change of the rotational 
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motion at the time of the collision and a short damped oscillation according to the excited 
inner degrees of freedom of the grain (distortion of beams) and their relaxation. After the 
collision the particles move on with constant translational and rotational velocities. For low 
impact rate we find a quite different behavior: First the particles collide as in the previous 
case, resulting in a certain rotational motion. Short time after the first collision, however, 
they collide a second time with different edges. Finally we find for the latter case of slow 
collisions a very small resulting angular velocity because the second collision causes angular 
moments in opposite directions, hence the resulting moment is small. 

Observing the motion of a moving particle colliding with a resting one (see [68]), i.e. the 
simplest possible contact of particles, one can remark that even for this system there is a 
complex behavior. For rigid spheres the velocities after a collision as a function of the impact 
velocities (eq. (1)) can be calculated analytically solving the viscoelastic equations of the 
spheres [8]. We guess that this will not be possible for the case of our particles. 

B. Outflow of a hopper 

The outflow of a hopper is of high interest not only because of the technological impor- 
tance of this process, but also because of the exciting phenomena as clogging and density 
waves which have been discovered recently (e.g. [69]). Density waves and clogging have 
been investigated using molecular dynamics in two and three dimensions in various papers, 
e.g. [23,70,71]. There is at least one other interesting and surprising phenomenon, detected 
recently by Evesque and Meftah [72,73]: They found that a hour glass "ticks" much slower 
if it is subjected to vertical vibrations y = A cos (cot). The effect seems to depend on the 
frequency / = u/2 7r and on the acceleration of the vibration T = Alo 2 . For frequencies 
between 40 Hz < f < 60 Hz surprisingly they observed for some accelerations that the flow 
almost stops. 

We investigated this effect using different particle models and we found that the effect 
could be reproduced with the new particle model only. The results will be described in 
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detail in a forthcoming paper [74]. 

Fig. (17) shows snapshots of the system with N = 400 complex particles. For an animated 
sequence we refer to our WWW-site [68]. The left figure shows the outflowing hopper, the 
flow varies irregularly and the grey scale codes for the particle velocities. When animating 
the figures one clearly observes density waves [68]. The right hand figure shows the same 
system a few moments later. The flow has stopped due to clogging. 

C. Granular flow in a rotating cylinder 

The flow of granular material in a rotating cylinder is one of the most popular problems 
in the field of granular materials. It has been investigated experimentally and theoretically 
by many authors using various techniques (see e.g. [75-78,27,43,50] and many others). The 
results shall not be discussed here. We applied the algorithm to this problem too. For the 
detailed description of the results see [79] . Here we only want to demonstrate the abilities 
of the method and to mention the main results briefly. 

Fig. 18 shows snapshots of the simulation of a slowly rotating cylinder. The grey scale 
codes for the velocity of the grains, black means high velocity, black codes for high velo- 
city. For low angular velocity of the driven cylinder one finds experimentally stick-slip flow 
(e.g. [75]), i.e. the material moves downwards not homogeneously but it forms avalanches. 
In the right figure one observes an avalanche on the top of the material indexed by light grey 
shadowed grains. The left figure shows the same system immediately before the avalanche. 
Avalanches are relatively seldom events. Most of the time the particles rest with respect 
to each other, i.e. they move only due to the rotation of the cylinder (left side of fig. 18). 
When increasing the angular velocity of the cylinder there is a relatively sharp transition 
between the stick-slip flow and the homogeneous regime. This transition could be found 
in the simulation using non-spherical grains [79]. Animated sequences for both regimes, 
stick-slip and continuous flow, can be accessed via World Wide Web [68]. 

Fig. 19 shows the dependence of the difference between the inclination of the material 
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and the angle of repose on the angular velocity of the cylinder. The line displays the behavior 
found by Rajchenbach 9 — G c ~ Q 2 [75]. The simulation data are in good agreement with 
this observation. 

Both results mentioned, the sharp transition between the flow regimes and the inclination 
as a function of the angular velocity of the cylinder, agree with experimental observations. 
We have not been able to find these effects using simple spheres or using particles composed 
of spheres [50,29] in molecular dynamics simulations. 

VII. CONCLUSION 

We presented a model for the simulation of the particles of a granular material. Each 
particle consists of triangles which are connected by deformable beams. There are no rest- 
rictions concerning the shape of the grains (convex or concave), nor the number or the shape 
of the triangles a grain consists of. To preserve calculation time it is favourable to chose 
arrangements of the triangles so that the beams are fixed at the center of mass points of the 
triangles. Otherwise one needs additional computation time caused by the Steiner law. 

When two triangles of different grains collide there acts an elastic restoring force pro- 
portionally to the compression of the particles. During the collision of the triangles the 
energy is preserved. A collision of a triangle with another one causes moments and forces 
acting on this triangle and hence a deformation of the beams which connect the triangles to 
neighboring ones. The deformation of the beams causes forces and moments acting on the 
neighboring triangles. 

Beams can be deformed in axial direction (elongation) and in shear direction and they 
can be bent. When a beam is deformed it dissipates energy proportionally to the deformation 
rate. 

In the present paper we described beams that recover completely while dissipating me- 
chanical energy when the deforming moments and forces vanish. Our model, however, is 
not restricted to this case. Interesting phenomena are plastic deformation and wear which 
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can be easily simulated by introducing thresholds for the beam forces and moments or for 
the deformations. When a force or a moment or a deformation, respectively, exceeds the 
threshold the beam breaks or deforms permanently. Similar simulations with other models 
have been reported in the literature [55-57]. 

The proposed algorithm has been implemented in FORTRAN. We have demonstrated 
the behavior of a system of grains applying it to three examples of granular assemblies. In a 
simple collision simulation of two grains we discussed the trajectories of the grains as well as 
the evolution of the kinetic and rotational energies depending on the impact rate. Sample 
results for the flow out of a hopper and the motion of granular material in a rotating cylinder 
have been presented. In the latter case we found quantitative agreement of the simulation 
with the experiment that could not be found so far, neither with spherical grains nor with 
non-spherical grains of other type. 
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APPENDIX A: SHEAR DEFORMATION CAN BE EXPRESSED BY BENDING 

The deformation of a beam in the two dimensional space has three degrees of freedom: the 



distance of the relative positions of the ends of the beams L = \]{xa — x b ) + (va — Vb) , 
and the angles 0^ and G#. These values are the parameters in eqs. (32) and (33). The 
shear parameter A in eqs. (22) is a fourth parameter, i.e. one of them must be redundant. 
In the following we will show that shear deformation can be expressed by bending, at least 
in our approach of a linear deformation-force relation. 

Suppose we have a deformed beam as drawn in Fig. 14 (bottom). Then applying eqs. (22) 
we find for the forces and moments according to shearing 
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6EI 
12EI 



M s a = M s b = 7 —^A (Ala) 



F l = -F'b = 77771 A (Alb) 



and according to bending (eqs. (32) and (33)) 



9 FT 

M>=-— (A2a) 
4EI 

M b B = -—Q* B (A2b) 
F b A = -F B = --^e* B . (A2c) 



The total forces and moments read 



6EI . 2EJ 



Mi = M% + M» = ^ A - — G^ (A3a) 

M B = M B + M b B = ^ A-^fe B (A3b) 

F* A = -F* = l -^A-^e B . (A3c) 
A 5 (L*) (L*) 5 V ; 

Turning the entire system by the angle G^ ~ tanG^ = — we obtain the new system as 
drawn in the top of fig. 14 with the new angles G^ = — ^ and G# = <d* B — ^. With eqs. (32) 
and (33) we find 

AEI 2EI 6EI 2EI ± . . 

M A = -—e A - — e B = ^ r A- — e* B (A4a) 
M^-^e.-^e^^A-^es, (A4b) 

F A = -F B = -™ { e A + e B ) = ^-^& B (A4c) 

Approximating L by V-7 2 + A 2 = L* evidently the results in eqs. (A4) coincide with 
eqs. (A3). Thus we do not need to consider shear deformation in our simulations since it is 
redundant here. 
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TABLES 



Basic algorithm 


portion of 


subroutine 


portion of 




computation 




computation 


Construction of the lists 


73.7 % 


Verlet list 


2.07 % 






Classifying interaction 


70.62 % 


Calculation of the 


7.87 % 


TYPE 1 


7.54 % 


triangle interaction 




TYPE 2 


0.24 % 






TYPE 3 


0.09 % 


Calculation of the 






13.33 % 


beam forces 








integration 


3.71 % 


predictor 


1.16 % 






corrector 


2.55 % 


all others 






1.4 % 



TABLE I. The performance analysis of the algorithm. The calculation has been done on a 
DEC-3000/700 workstation (explanation in the text). 



31 



FIGURES 





FIG. 1. Examples of grains composed of different numbers of triangles. The model is not 
restricted to convex grains. 




FIG. 2. The "standard" type of a collision between two triangles AiBiCi and A k B k C k . The 



forces ±r ifc are directed perpendicular to the intersection line SiS 2 - Moments act with respect to 



the center of mass points G { and Gj. The absolute value of the interaction force 
to the shadowed intersection area A (SiS 2 A k ). 



is proportional 
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A k A k 

FIG. 3. The first type of collisions: There are two intersection points Si and S 2 . Both lie at 
the same edge of one of the triangles. 




FIG. 4. The second type of collisions: There are two intersection points Si and S 2 lying on 
different edges for both triangles. 
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FIG. 5. The third type of collisions: There are four intersection points S±, S 2 , S 3 and S 4 . Each 
pair of points lies on one edge of either the first or the second triangle. 
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FIG. 7. The fifth type of collisions: There are six intersection points S 



lv . . ,0 6 . 




FIG. 8. An infinitesimally deformed beam with radius of curvature p. The fibres below or above 
the neutral fibre SS are in tension or in compression. For a beam with quadratic cross section 6 2 



the deformation leads to the resulting moment M 
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FIG. 9. The shape and the location of a beam at rest. Its length is assumed to be L, its ends 
A and B lie on the x-axis. When the beam is not deformed there do not act any forces F A , F B 
or moments M A , M B . In the following figures 8-11 and in fig. 14 the a;-axis is drawn as a thin 
horizontal line. 




FIG. 10. When the beam in fig. 9 undergoes a shear deformation A there act the forces F A 
and F B and the moments M A and M B . The forces and moments are given in eqs. (22). 
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FIG. 11. A bent beam. The forces F A and F B and the moments M A and M B can be calculated 
by superposing the forces and moments of two single sided bent beams (fig. 12), one of them bent 
at the side A, the other one at the side B. The results are given in eqs. (32) and (33). 




FIG. 12. A single sided bent beam. The forces and moments result from the superposition of 
the deformation drawn in fig. 13 and the restoring moment M B (fig. 12) which assures the angle 
@b to be zero. The equations for the forces and moments are given in eqs. (27) and (29). 




L 

FIG. 13. If the ends of the beam are not fixed but freely moving, the applied moment M causes 
bending due to the angles @ A and @ B given in eqs. (23). 



36 




FIG. 14. The deformation of a beam in shear direction A (upper part) can be represented by 
bending deformation (lower part). The appendix contains the proof that one need not to consider 
shear and bending deformation in the linear approximation which is justified for small deformations. 
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FIG. 15. The collision of a horizontally from left to right moving grain (black drawn) with a 
resting one (grey) for different impact rates v re i = 10 cm/sec (top), v re i = 50 cm/sec (middle), and 
v re i = 100 cm/sec (bottom). The stroboscopic snapshots of the grains are take at equidistant time 
intervals. An animated sequence of the collisions can be accessed via World Wide Web [67]. The 
right hand figures display the kinetic, rotational and total energies of the particle collisions over 
time. 
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FIG. 16. The angular velocities of the particles over the time. For lower velocities the particles 



collide two times, therefore those collisions results in a small angular velocity of the grains. 
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FIG. 17. Snapshots of the simulation of an outflowing hopper. The left figure shows regular 
outflow, the right figure show a clogging situation. An animated sequence is available via World 
Wide Web [67]. 




FIG. 18. Snapshots of the simulation with N = 500 complex particles with P G [0.1 cm; 0.2 cm] 
in a rotating cylinder of diameter D = 4 cm. The angular velocity is Q = 0.1 sec -1 . The wall 
particles have not been drawn. The right snapshot shows an avalanche, the left one has been taken 
short time before. The grey scale codes for the particle velocity. 
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FIG. 19. The inclination of the material surface over the angular velocity 0. The dotted line 
displays the function which has been measured experimentally by Rajchenbach [73]. 
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